DALS027-批次效应03-校正批次效应

发表于 | 分类于
| 字数统计: 2,986 | 阅读时长 ≈ 13

前言

这一部分是《Data Analysis for the life sciences》的第10章批次效应的第3小节，这一部分的主要内容涉及批次效应(Batch Effects)的校正，有关这一部分Rmarkdown文档参见作者的Github。

使用情境

在这一部分里，我们要了解如何对批次效应进行校正。

数据案例

为了说明我们可以使用统计方法来校正批次效应，我们将创建一个数据示例，其中感兴趣的结果与批次之间有某种程度的混杂，但不完全混杂。为了辅助说明和评估我们所示的方法，我们还将选择一个结果用于说明我们期望对哪些基因应该是差异表达的。也就是说，我们想让性别作为我们感兴趣的变量，并且期望 Y 染色体上的基因表达有所差异。我们或许能看到，来源于 X 染色体上的基因是差异表达的，因为有些基因在 X 染色体上并没有失活，而是正常表达，这些数据都包含在以下的数据集中，如下所示：

1
2
3
##available from course github repository
library(GSE5859Subset)
data(GSE5859Subset)

我们可以看到，性别与月份有相关性：

1
2
month <- format(sampleInfo$date,"%m")
table(sampleInfo$group, month)

结果如下所示：

1
2
3
4
5
6
> month <- format(sampleInfo$date,"%m")
> table(sampleInfo$group, month)
   month
    06 10
  0  9  3
  1  3  9

为了说明混杂现在，我们会挑选一些基因绘制热图。现在我们挑选这些基因：(1)所有Y染色体上的基因；(2)某些与批次效应有效的基因；(3)一些随机选择的基因。下面是热图，其中红色表示高表达基因，蓝色表示低表达基因，黄色表示中间表达基因。每列是一个样本，每行就是随机选择的基因，代码如下所示：

image_of_subset, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
library(rafalib)
library(RColorBrewer)
library(genefilter)
batch <- factor(format(sampleInfo$date,"%m"))
chr <- geneAnnotation$CHR
tt<-rowttests(geneExpression,batch)
ind1 <- which(chr=="chrY") ##real differences
ind2 <- setdiff(c(order(tt$dm)[1:25],order(-tt$dm)[1:25]),ind1)
set.seed(1)
ind0 <- setdiff(sample(seq(along=tt$dm),50),c(ind2,ind1))
geneindex<-c(ind2,ind0,ind1)
mat<-geneExpression[geneindex,]
mat <- mat -rowMeans(mat)
icolors <- colorRampPalette(rev(brewer.pal(11,"RdYlBu")))(100)
mypar(1,1)
image(1:24,1:nrow(mat), t(mat),xaxt="n",yaxt="n",col=icolors,xlab="",ylab="")
axis(3,1:24,rep(c("F","M"),each=12),cex.axis=0.5)
axis(1,1:24,month,cex.axis=0.5)

在上图中，前12列是女性（用1表示），后12列是男性（用0表示）。我们可以看到，一些Y染色体上的基因趋向顶部，其中在女性数据中这些是蓝色，在男性数据中这是红色。我们也可以看到一些基因与月份相关，这些基因集中趋向于底部。一些基因在六月份低表达（June，6），一些基因在10月份高表达（October，10），还有一些基因正好相反。月份效应并没有性别效应明显，但是这种效应确实存在。

在接下来的内容里，我们会模拟实际分析中的一些典型分析方法。我们会假设一种情况，在这种情况里，我们假设不知道哪些基因在男女性别有差异，然后我们要找出这些差异基因，通过比较我们所期望的相关性，来评估一下这些方法的优劣。这里需要注意一下，我们虽然只绘制出其中一部分基因，但是我们所分析的基因数目却是8793个。

评估图与总结

为了说明我们所用方法的优劣，我们将假设常染色体上的基因（这些基因不在染色体X和染色体Y上）有可能是假阳性，位于Y染色体上的基因是真阳性。染色体X上的基因可以是假阳性，可以是真阳性。上述的这些设计就在让我们评估特异性(specificity)与灵敏性(sensitivity)。由于在实际分析中，我们很少能知道“真相”，因此我们在这里做的评估是很在实际分析中实现的。因此模拟常常用于评估目的：我们知道真相，因为分析所用的数据是我们自己构建的。但是，模拟也会有无法捕获真实实验数据细微差异的风险。相比之下，这个数据集就是一个实验性质的数据集。

在接下来的部分中，我们将使用p值的直方图来评估一下批次效应校正程序的特异性(specificity)，特异性表示低假阳性率。因为常染色体应该不会出现由于差异导致的差异表达，我们应该会看到一个平坦的p值直方图。为了评估灵敏性(sensitivity,低假阳性)，我们会报告那些位于X染色体和Y染色体上我们拒绝零假设时，所计算出的基因的数目。我们还会生成一个火山图，用一条水平虚线将那些有显著性差异的基因与无显著性差异的基因区分开来，用于突出显示X染色体和Y染色体的基因。以下是一个应用原始t检验得到的结果，以及q值小于0.1的基因，如下所示：

pvalue_hist_and_volcano_plots, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(qvalue)
res <- rowttests(geneExpression,as.factor( sampleInfo$group ))
mypar(1,2)
hist(res$p.value[which(!chr%in%c("chrX","chrY") )],main="",ylim=c(0,1300))
plot(res$dm,-log10(res$p.value))
points(res$dm[which(chr=="chrX")],-log10(res$p.value[which(chr=="chrX")]),col=1,pch=16)
points(res$dm[which(chr=="chrY")],-log10(res$p.value[which(chr=="chrY")]),col=2,pch=16,
       xlab="Effect size",ylab="-log10(p-value)")
legend("bottomright",c("chrX","chrY"),col=1:2,pch=16)
qvals <- qvalue(res$p.value)$qvalue
index <- which(qvals<0.1)
abline(h=-log10(max(res$p.value[index])))
cat("Total genes with q-value < 0.1: ",length(index),"\n",
    "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
    "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")

很明显，我们注意到直方图并不平坦的。相反，低p值的条状过高（原文：low p-values are over-represented）。此外，最终列表上超过一半的基因是常染色体。我们现在介绍两个统计学方法来尝试一下解决这个问题。

使用线性模型来校正批次效应

我们已经观察到了处理时间(processing date)对基因表达有着某种效应。我们因此会通过将此效应添加到一个模型中来对其进行校正(adjust)。当我们对两组数据进行t检验时，它就相当于对基因 $j$ 进行以下线性模型拟合：

其中，研究对象 $i$ 如果是女性，那么 $x_{i}=1$, 如果是男性则为0，$\beta_{j}$ 表示基因 $j$ 的差异，$\varepsilon_{ij}$ 表示组内变异。现在我们的问题是什么？

在线性模型那一章节中，我们假设错误项(error term)是独立的。我们知道，对于所有的基因来说，这是不可能的，因为我们知道10月份数据中的错误项比6月份中的错误项更相似。我们可以在线性模型中添加一项，用于校正这种错误，如下所示：

在这个公式中，如果样本 $i$ 是在10月份处理的，那么 $z_{i}=1$ ，否则为0，另外， $\gamma_{j}$ 是基因$j$的月份效应。这就是 一个线性模型如何比t检验具有更大灵活性的一个例子。

我们构建一个含有批次的模型矩阵，如下所示：

1
2
sex <- sampleInfo$group
X <- model.matrix(~sex+batch)

现在我们对每个基因拟合一个模型。例如，我们要注意在原始模型和那些含有批次校正模型之间的差异，如下所示：

1
2
3
4
5
6
7
8
j <- 7635
y <- geneExpression[j,]
X0 <- model.matrix(~sex) 
fit <- lm(y~X0-1)
summary(fit)$coef
X <- model.matrix(~sex+batch)
fit <- lm(y~X)
summary(fit)$coef

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
> summary(fit)$coef
                Estimate Std. Error   t value     Pr(>|t|)
X0(Intercept)  6.9555747  0.2166035 32.112008 5.611901e-20
X0sex         -0.6556865  0.3063237 -2.140502 4.365102e-02
> X <- model.matrix(~sex+batch)
> fit <- lm(y~X)
> summary(fit)$coef
               Estimate Std. Error    t value     Pr(>|t|)
(Intercept)  7.26329968  0.1605560 45.2384140 2.036006e-22
Xsex        -0.04023663  0.2427379 -0.1657616 8.699300e-01
Xbatch10    -1.23089977  0.2427379 -5.0709009 5.070727e-05

我们然后将这个新模型用于每个基因，例如我们可以使用sapply()函数来恢复估计的系数和p值，如下所示：

pvalue_hist_and_volcano_plots2, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
res <- t( sapply(1:nrow(geneExpression),function(j){
  y <- geneExpression[j,]
  fit <- lm(y~X-1)
  summary(fit)$coef[2,c(1,4)]
} ) )
##turn into data.frame so we can use the same code for plots as above
res <- data.frame(res)
names(res) <- c("dm","p.value")
mypar(1,2)
hist(res$p.value[which(!chr%in%c("chrX","chrY") )],main="",ylim=c(0,1300))
plot(res$dm,-log10(res$p.value))
points(res$dm[which(chr=="chrX")],-log10(res$p.value[which(chr=="chrX")]),col=1,pch=16)
points(res$dm[which(chr=="chrY")],-log10(res$p.value[which(chr=="chrY")]),col=2,pch=16,
       xlab="Effect size",ylab="-log10(p-value)")
legend("bottomright",c("chrX","chrY"),col=1:2,pch=16)
qvals <- qvalue(res$p.value)$qvalue
index <- which(qvals<0.1)
abline(h=-log10(max(res$p.value[index])))
cat("Total genes with q-value < 0.1: ",length(index),"\n",
    "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
    "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")

结果如下所示：

从图上我们可以看到，特异性（假阳性较少）有了很大的提高，同时灵敏性也没有过多损失（我们仍然能发现许多Y染色体上的基因）。但是，我们仍然能够从直方图中看到一些偏差，在后面的部分里，我们可以看到月份并不能完美地解释批次效应，应该有更好的估计方法。

注意计算效率

在上面代码中，当我们重复计算 $(X^\top X)^{-1}$，以及将其应用于每个基因时，设计矩阵并没发生改变。相反，我们可以在一次矩阵代数计算中实现该计算，然后通过将 $(X^\top X)^{-1}X^\top Y$ 与表示所有基因的 $Y$ 列相乘来获得所有的 $\beta$ 。limma包的设计就是基于这个思路（利用了QR分解）。我们来看一下这样的计算有多快：

1
2
3
library(limma)
X <- model.matrix(~sex+batch)
fit <- lmFit(geneExpression,X)

每个基因的估计回归系数如下所示：

1
2
3
dim( fit$coef)
# > dim( fit$coef)
# [1] 8793    3

我们对每个基因都有一个估计。如果要获得其中一个p值，我们就必须要构建以下的比值(ratio)：

We have one estimate for each gene. To obtain p-values for one of these, we have to construct the ratios:

1
2
3
4
k <- 2 ##second coef
ses <- fit$stdev.unscaled[,k]*fit$sigma
ttest <- fit$coef[,k]/ses
pvals <- 2*pt(-abs(ttest),fit$df)

Combat

Combat 是一种基于线性模型的算法，常用于校正批次效应。这种方法能够对估计值拟合一个分层模型，并且移除行的特定批次效应。Combat使用了一个模块块的方法(a modular approach)。第一步就是删除那些被认为是批次效应的内容，如下所示：

message
1
2
3
library(sva) #available from Bioconductor
mod <- model.matrix(~sex)
cleandat <- ComBat(geneExpression,batch,mod)

然后，将剩下的结果用于拟合我们感兴趣的变量：

1
res<-genefilter::rowttests(cleandat,factor(sex))

在这种情况下，与我们通过拟合简单线性模型得到的结果相比，通过上面计算得到的结果缺乏特异性（原文：the results are less specific than what we obtain by fitting the simple linear model）：

pvalue_hist_and_volcano_plots3, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
mypar(1,2)
hist(res$p.value[which(!chr%in%c("chrX","chrY") )],main="",ylim=c(0,1300))
plot(res$dm,-log10(res$p.value))
points(res$dm[which(chr=="chrX")],-log10(res$p.value[which(chr=="chrX")]),col=1,pch=16)
points(res$dm[which(chr=="chrY")],-log10(res$p.value[which(chr=="chrY")]),col=2,pch=16,
       xlab="Effect size",ylab="-log10(p-value)")
legend("bottomright",c("chrX","chrY"),col=1:2,pch=16)
qvals <- qvalue(res$p.value)$qvalue
index <- which(qvals<0.1)
abline(h=-log10(max(res$p.value[index])))
cat("Total genes with q-value < 0.1: ",length(index),"\n",
    "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
    "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")

练习

P442